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SUMMARY 

This  paper  presents  two  finite  element  approaches  for  the  three  dimensional  elastic- 
plastic  analysis  of  a  surface  crack.  The  distributions  of  the  stress  and  strain  fields  around 
the  crack  are  discussed  in  detail  and  it  is  shown  that  both  approaches  yield  similar  values 
for  the  strain  energy  distribution  around  the  crack. 
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Vector  of  nodal  displacements 

System  of  constraint  equations 

Path  independent  line  integral  around  crack  tip 

Plastic  stress  intensity  factor 

Plastic  strain  intensity  factor 

Penalty  element  stiffness  matrix 

Matrix  of  coefficients  in  system  of  constraint  equations 

Exponent  in  power  law  hardening  equation 
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Displacement  component 

Strain  energy  density 

Local  isoparametric  coordinates 

Cartesian  co-ordinate  system 

Constant  in  power  law  hardening  equation 

Strain  tensor 

Reference  value  of  strain 

Stress  tensor 

Reference  value  of  stress 

Global  potential  energy  functional 


1.  INTRODUCTION 


In  recent  years  a  great  deal  of  work  has  been  undertaken  in  order  to  characterize  the  elastic 
and  plastic  stress  strain  fields  around  a  crack  which  is  either  in  plane  strain  or  plane  stress  [1-4]. 
A  number  of  criteria  for  assessing  severity  have  been  proposed;  of  these  the  J  integral  approach 
[3]  and  the  strain  energy  density  factor  (viz.,  S)  approach  [S]  are  widely  used  and  it  has  been 
shown  that  S  and  J  are  linearly  proportional. 

The  purpose  of  the  present  paper  is  to  develop  a  methodology  for  three  dimensional  elastic- 
plastic  finite  element  analysis  and  to  investigate  whether  the  functional  form  of  the  two 
dimension  solution  is  valid  for  a  three  dimensional  crack. 


2.  TWO  DIMENSIONAL  ANALYSIS 

Let  us  consider  a  material  for  which  the  relationship  between  stress  and  strain  in  simple 
tension,  is  given  by: 


«/«0  =  a(o/oo)”  (2.1) 

where  a  is  a  dimensionless  constant  and  er0  and  <„  are  reference  values  of  the  stress  and  strain. 
In  plane  strain  the  stress,  strain  and  displacement  fields  of  the  dominant  singularity  at  the  crack 
tip  have  the  form: 


[oij,  ®e]  =  <7oAf„r')'”+1[Sfj(0),  8^8) ] 

(22) 

<11  =  x«oA>nnUiij(0) 

(2.3) 

ui  =  atoAer1/n+li}|(0) 

(2.4) 

where  r  is  the  distance  from  the  crack  tip  and  8  is  the  angle  measured  from  directly  ahead  of  the 
crack  tip.  The  dimensionless  functions  an,  oe,  cij  and  lit  depend  on  the  strain  hardening  exponent 
n  and  have  been  given  in  [I],  The  plastic  stress  and  strain  intensity  factors  K„  and  K ,  (=  A„") 
are  related  to  the  J  integral  by: 

J  —  auQfoK^Kel 


=  xa0eoA„n"/ 


where  /  is  only  a  function  of  n  and  is  tabulated  in  [I],  The  strain  energy  density  It' in  the  direction 
8  —  0,  i.e.  ahead  of  the  crack,  defined  by  : 


W  -  I  xtjrfcij 

is  related  to  K ,  and  K„  by 


(2.6) 


(2.7) 


I 


so  that  the  strain  energy  density  function  S, 
equation: 

s  =  Jn  . 

/(«+ 1) 


where  S  =  rW,  can  be  directly  related  to  J  by  the 


(2.8) 


In  deriving  equations  (2.2),  (2.3)  and  (2.4)  references  (1-4]  used  several  simplifying  assump¬ 
tions.  Indeed,  these  equations  are  only  valid  for  monotonic  loading.  Nevertheless,  this  repre¬ 
sentation  is  commonly  used. 

If  the  stress-strain  relationship  is  assumed  to  be  piecewise-linear,  rather  than  as  given  in 
equation  (2.1),  then  the  stresses  exhibit  an  r1  variation  at  the  crack  tip  [1],  rather  than  the 
r-i  n+i  variation  mentioned  previously.  However,  predictions  for  the  strain  energy  density  M'', 
and  hence  S  and  J,  at  the  crack  tip  are  relatively  insensitive  to  the  assumed  stress-strain 
relationship. 

For  the  special  case  when  n  —  1  and  the  material  is  purely  elastic,  the  above  equations 
reduce  to  the  well  known  equations  for  two-dimensional  linear  elastic  fracture  mechanics.  In 
linear  elastic  fracture  mechanics  theory  the  functional  form  of  the  stress  field  around  a  three 
dimensional  crack  is  the  same  as  that  around  a  crack  in  two  dimensional  plane  strain.  This 
raises  the  possibility  that  in  three  dimensional  elastic-plastic  fracture  mechanics,  the  functional 
form  of  the  stress  and  strain  fields  may  be  similar  to  that  given  in  equations  (2.2),  (2.3)  and  (2.4). 


3.  A  PENALTY  FOR  ELASTIC-PLASTIC  FRACTURE 

Let  us  now  assume  that  the  displacements  around  a  three  dimensional  crack  have  the 
functional  form  given  in  equation  (2.4).  This  assumption  can  be  built  into  the  standard  finite 
element  method  in  a  variety  of  ways,  the  simplest  of  which  is  the  penalty  approach  [6],  which 
we  will  now  term  method  A.  In  this  approach  we  minimize  the  constrained  functional  II*  (see 
Appendix)  where 

3 

11*  --  I  I  4-  P 

i  =  I 

and  where  II  is  the  complementary  potential  energy.  P  is  the  penally  number,  W'(.v,  y,  z)  is 
a  positive  weighting  function,  ut  are  the  components  of  the  displacement  field  around  the  crack 
for  a  standard  isoparametric  element,  whilst  ut  are  the  expressions  for  the  near  tip  displacement 
field  as  given  in  (2.4). 

When  P  becomes  large  the  constraint  u i  —  u\  is  enforced.  In  this  work,  fifteen-noded 
isoparametric  wedge  elements  were  used  around  the  crack  tip  and  different  weighting  functions 
evaluated.  However,  based  on  physical  grounds,  the  approach  adopted  was  to  use  reduced 
integration,  as  required  in  [6],  and  set 

W  =  I  at  the  nodes  and  at  the  near  tip  Gauss  points 

=  0  at  all  other  points. 

In  the  next  section,  the  validity  of  this  approach  will  be  evaluated  by  comparing,  amongst  other 
things,  the  values  for  the  strain  energy  density  function  S  computed  from  the  elements  in  front 
of  the  crack  tip  elements,  with  the  value  of  5  computed  using  the  relationship 

5  =  <xo0f0  A„n  1 1 
n+  I 

,  «ao«o  (3.2) 

n+  I 

where  K,  is  evaluated  from  the  opening  of  the  crack,  i.e.  using  equation  (2.4). 


rr  • 

v  v  a 


W(x,  y,  zXui-ui  )zdv 


(3.1) 


3.1  Alternative  Computational  Methods 


A  variety  of  methods  have  been  developed  for  the  finite  element  analysis  of  two  dimensional 
elastic-plastic  fracture  problems.  The  best  of  these  approaches  are  reviewed  in  [7,  8]  and 
involve  using  either: 

(i)  crack  tip  elements  with  the  midside  nodes  moved  to  the  quarter  points  and  with  the 
stiffness  matrix  computed  using  reduced  integration; 

(ii)  crack  tip  elements  with  the  midside  nodes  moved  to  the  quarter  points  and  all  nodes 
at  the  crack  tip  allowed  to  have  different  degrees  of  freedom; 

(iii)  very  small  elements  at  the  crack  tip  with  each  element  having  a  Young's  modulus  of 
unity  (this  procedure  is  known  as  Unimod). 

Of  these  approaches,  the  first,  which  we  will  now  term  method  B,  is  the  simplest  and  yields  very 
accurate  results  [8].  The  other  methods  are  easy  to  implement  for  two  dimensional  problems  but 
are  very  tedious  to  implement  for  three  dimensional  problems.  The  second  approach  also 
generates  an  r  1  strain  singularity  which  is  only  true  for  an  elastic  perfectly  plastic  material. 


4.  NUMERICAL  INVESTIGATION 


Let  us  now  consider  the  applicability  of  the  approaches  described  above  for  the  solution 
of  a  three  dimensional  elastic-plastic  fracture  problem.  Here,  we  analyse  the  distribution  of 
the  strain  energy  density  W,  and  .S  around  a  part  circular  crack  in  an  axially  loaded  specimen. 
The  specimen  is  shown  in  Figure  I  and  was  chosen  to  coincide  with  the  specimens  tested  in 
reference  [9],  The  material  was  taken  to  obey  the  stress-strain  law. 


€ 


(4.1) 


with 


and 


e0  =  2- 175 x  10  3 


<t0  —  443  -4  MPa. 


Because  of  the  symmetrical  nature  of  the  problem  only  1/4  of  the  structure  was  modelled, 
see  Figure  2.  The  resultant  finite  element  mesh  consisted  of  12  of  the  15-noded  isoparametric 
wedge  elements  and  249  of  the  20-noded  isoparametric  brick  elements.  The  length  scale  for  the 
crack  tip  elements  was  1/50  of  the  tadius  of  curvarure  of  the  crack  front.  Two  different  analyses 
were  conducted.  Firstly,  in  what  we  term  method  A.  the  crack  tip  elements  used  were  the  penalty 
elements  as  formulated  in  Section  3.  Secondly,  in  method  B,  the  crack  tip  elements  had  their 
midside  nodes  moved  to  the  quarter  point  positions.  The  stiffness  matrices  for  all  of  the  elements 
were  evaluated  in  double  precision  using  reduced  integration.  The  maximum  applied  stress 
considered  was  177  MPa  and  at  each  load  level  the  load  steps  were  varied  to  ensure  that  conver¬ 
gence  had  been  obtained.  Solutions  were  performed  using  the  PAFEC  suite  of  finite  element 
programs,  on  the  ARL  VAX  1 1/780  computer. 

In  order  to  compare  these  two  approaches,  the  values  of  S(=  rW)  around  the  crack  front, 
are  given  in  Table  I  at  three  different  load  levels.  Here  S  is  computed  from  the  elements  directly 
in  front  of  the  crack. 


TABLE  I 


Strain  Energy  Density  Function  S  Around  Crack  Front 


Applied  Stress 

a 

MPa 

Method 

Strain  Energy  Density  Function 

S 

(N/m) 

Position  Around  Crack 

1 

2 

3 

4 

5 

6 

88  6 

A 

230 

221 

241 

253 

333 

380 

B 

234 

244 

258 

283 

370 

459 

132  9 

A 

879 

884 

903 

935 

1200 

1362 

B 

843 

884 

899 

945 

1188 

1581 

177-2 

A 

2235 

2100 

2235 

2287 

2875 

3257 

B 

2088 

2166 

2126 

2*50 

2659 

3874 

With  the  exception  of  point  P6  which  is  a  near  surface  point,  both  approaches  give  very 
similar  estimates  for  5.  Indeed  as  the  load  increases  the  discrepancies  in  these  two  approaches 
reduces.  This  is  to  be  expected  since  at  the  lower  loads  not  all  of  the  Gauss  points  in  the  crack 
tip  elements  had  gone  plastic. 

To  assess  the  validity  of  equations  (2.2),  (2.3)  and  (2.4),  which  are  commonly  referred  to 
as  the  HRR  equations,  S  was  also  computed  indirectly,  using  equation  (3.2),  from  the  opening 
of  the  crack.  Indeed  these  values,  divided  by  the  value  of  S  at  point  PI.  are  shown  in  Table  2. 


TABLE  2 


Factorised  Strain  Energy  Density  Function  S  Around  Crack  Front 


Applied 

Stress 

rj 

MPa 

Method 

Factorised  Strain  Energy  Density  Function  S/S pI 

Position  Around  Crack 

1 

2 

3 

4 

6 

88-6 

A  —direct  approach 

1  00 

0-96 

1  05 

110 

1  45 

1  65 

A  —indirect  approach 

1  -00 

1  -00 

1  01 

1  05 

1  06 

1  39 

B — direct  approach 

1-00 

1  -09 

1-15 

1-26 

1  65 

2  05 

B — indirect  approach 

I  00 

0-99 

1  01 

1  03 

1  -23 

I  3) 

132  9 

A — direct  approach 

t  -00 

I -04 

1  07 

1  - 12 

1  41 

1  88 

A — indirect  approach 

1  00 

0  99 

0-98 

0  99 

1  01 

1  32 

B  — direct  approach 

1  00 

0-95 

1  -03 

1  07 

1  -37 

1  55 

B— indirect  approach 

1-00 

0-99 

0-99 

1  02 

1  02 

1  -29 

177-2 

A— direct  approach 

1-00 

0-94 

1  -00 

1-02 

1  29 

1  45 

A— indirect  approach 

1  00 

0-98 

0  97 

0  96 

0  96 

I  -27 

B — direct  approach 

1-00 

1  -04 

1  02 

1  03 

1  -27 

I  85 

B — indirect  approach 

1  -00 

1  01 

0  97 

1  02 

1-15 

1  29 

In  general,  the  variation  around  the  crack  is  similar  with  the  maximum  discrepancy  at  the 
lower  toad  levels.  As  previously  mentioned  not  all  of  the  Causs  points  in  the  crack  tip  elements 
were  plastic  at  the  lower  loads.  Since  the  stress-strain  law  used  in  the  numerical  analysis  was  a 
piece  wise  linear  approximation  to  equation  (4.1),  it  appears  that  equations  (2.2),  (2.3)  and  (2.4) 
represent  a  reasonable  first  approximation.  Indeed  assuming  failure  to  occur  when  S  =  Sc, 
where  Sc  is  the  critical  value  of  5  for  the  material,  the  maximum  discrepancy  in  the  computed 
values  of  S  corresponds  to  only  a  12%  difference  in  failure  load. 

The  computed  values  of  5  depend  more  on  the  way  in  which  it  is  evaluated,  i.e.  by  the  direct 
or  indirect  approach,  than  on  the  method  by  which  the  solution  was  obtained,  i.e.  method  A 
or  B.  The  simplest  method  is  to  use  isoparametric  elements  around  the  crack  with  their  mid¬ 
side  nodes  moved  to  the  quarter  points  and  evaluate  W  and  S  directly  from  the  elements  in  front 
of  the  crack  tip  elements.  Reduced  integration  is  recommended  when  formulating  the  stiffness 
for  the  elements  [8]. 


5.  CONCLUSION 

This  paper  has  shown  that  the  two  dimensional  elastic-plastic  f racture  equations  represent 
a  first  approximation  to  the  stress  and  strain  fields  around  a  three  dimensional  crack.  As  in 
two  dimensional  elastic-plastic  fracture  mechanics,  the  strain  energy  density  function  is  relatively 
insensitive  to  the  method  by  which  it  is  computed. 
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Footnote.  It  should  be  noted  that  the  failure  load  is  proportional  to  (S),'n+1. 
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APPENDIX 


A  Penalty  Eeimeat  for  Plastic  Fracture 

In  order  to  enforce  the  functional  form  of  displacements  around  the  crack  tip  we  use  the 
constrained  variational  approach  detailed  in  reference  [6].  Consider  a  typical  15  noded  isopara¬ 
metric  element  at  the  crack  tip,  as  shown  in  Figure  3.  We  wish  to  enforce  equation  (2.4)  for  the 
nodes  along  the  six  radial  lines.  For  simplicity  we  first  investigate  the  displacement  field  in  one 
dimension  only,  as  shown  in  Figure  4.  In  this  simplistic  approach  the  element  has  three  nodes 
along  its  length  /,  with  node  1  at  the  crack  tip.  We  require  displacements  to  be  in  accordance  with 
equation  (2.4),  so  that 


u'  =  ao+airl"t*l+air>  (Al) 

with  ao,  ai  and  a-i  as  arbitrary  constants  which  are  to  be  determined.  Application  of  the  boundary 
conditions;  (i)  u  =  ui  at  r  =  0,  (ii)  u  =  u2  at  r  =  1/2  and  (iii)  u  =  us  at  r  =  /,  to  equation  (Al) 
yields 


U  =  f\U\  +fiU2  +fsU3 


(A2) 


where 


. (;)'-(;) . ] 

tW"] 


f  _  _  V2 
J2,  |  _2<1-1/n+1 


/3- 


I 
I 


1M2-"- . -Hn-)] 


+  i 

(A3) 


However,  in  the  local  isoparametric  co-ordinate  system  X,  the  displacement  field  is  given  by 

u  =  -^(l-A')m  +  (l-A'8)(/2+^(l-A)«s  (A4) 

We  require  the  displacement  field  as  given  by  equation  (A2)  and  (A4)  to  be  identical.  Thus, 
the  constraint  condition  to  be  satisfied  is 

C(u)  =  u—u'  =  0.  (A5) 

Substituting  for  u  and  u  from  equations  (A2)  and  (A4)  respectively  into  equation  <A5) 
gives 

C(u)  =  f-iui-t-  LiUi  +  L3U3  (A6) 

where 

it 


6 


I 


and  fi,  f?  and  fz  are  as  in  equation  (A2). 

This  constraint  must  be  applied  to  all  radial  lines  in  the  penalty  element,  and  to  all  penalty 
elements  around  the  crack  front.  The  individual  constraint  equations  (i.e.  equation  (A6))  are 
combined  in  matrix  form  to  give 

C(u)  =  Uu)  (A7) 

where  L  is  the  matrix  of  constraint  coefficients,  and  u  is  the  vector  of  nodal  degrees  of  freedom. 

We  now  build  the  required  constraint  condition  C(u)  =  0,  of  equation  (A7)  into  the  standard 
finite  element  method  by  minimising  the  constrained  complementary  potential  energy  functional 
II*,  see  reference  [6], 


i 


w(x.y,z)C'T(u)C(u)dy 


(A8) 


which  appears  in  modified  form  in  Section  3.  Here  n  is  the  complementary  potential  energy, 
P  is  the  penalty  number,  and  w  is  a  positive  function.  As  the  penalty  number  P  is  increased 
the  constraint  condition  C(«)  =  0  is  enforced.  Substituting  for  C (u)  from  equation  (A7)  and 
differentiating  with  respect  to  the  nodal  degrees  of  freedom,  yields  the  penalty  stiffness  matrix 


2  P 


w{x,y,z)ULdV. 


(A9) 
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FIG.  1(a)  DIMENSIONS  OF  TEST  SPECIMEN 


FIG.  1(b)  DETAILS  OF  SPARK  MACHINED  CRACK  AT  SECTION  A-A. 


FIG.  4  CO-ORDINATES  OF  CRACK  TIP  ELEMENT  NODES 
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